knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE)

library(tidyverse)
library(skimr)
library(lubridate)
library(stringr)
library(nortest)
library(moments)
library(rstatix)
library(effectsize)
library(igraph)
library(tidygraph)
library(ggraph)
library(ggiraph)
library(broom)
library(robustbase)   
library(lmtest)       
library(performance)
library(see)
library(modelsummary)


# Formázott p-érték függvény
format_p <- function(p) {
  ifelse(p < 0.001, "< 0.001", as.character(round(p, 3)))}

Bevezetés

A gyógyszerek forgalomba hozatali engedélyezése összetett, több lépcsős szabályozási folyamat, amelyet az engedélyezést követően is folyamatos felügyelet és a dokumentáció ismtételt aktualizálása kísér. A forgalomba hozatali engedélyhez tartozó dokumentáció módosításai, az úgynevezett revíziók, a gyógyszer teljes életciklusa során tükrözhetik többek között a biztonsági jelzések, a klinikai bizonyítékok bővülésének, illetve a szabályozói eljárások sajátosságainak hatását. Emellett az engedélyezések mennyisége és összetétele (például terápiás területek szerint) időben is változhat, ami az gyógyszerfejlesztési trendek és a gyógyszercégek fókuszának változását mutathatják.

Az Európai Gyógyszerügynökség (European Medicines Agency, EMA) központi szerepet játszik a gyógyszerek engedélyezésében az Európai Unióban. Jelen elemzés célja az EMA gyógyszerengedélyezési adatbázisának részletes vizsgálata, amely négy fő kutatási kérdés köré szerveződik: 1. Az engedélyezések időbeli trendjeinek feltárása 1995 és 2023 között. 2. Az ATC-csoportok (Anatomical Therapeutic Chemical) közötti különbségek vizsgálata kulcsváltozók mentén 3. A revíziók számát befolyásoló tényezők azonosítása többváltozós modellezéssel 4. A központi idegrendszeri (ATC: N csoport) készítményeknél a terápiás területek és hatóanyagok kapcsolati szerkezetének hálózati szemléletű bemutatása.

Az elemzéshez az EMA nyilvánosan elérhető adatbázisát használtam, amely a TidyTuesday projekt keretében került publikálásra (2023. március 14.).

Adatok előkészítése

Az elemzés reprodukálható módon, webes forrásból betöltött adatbázison alapul, közvetlenül a TidyTuesday GitHub tárhelyéről került beolvasásra.

# Adat beolvasása a GitHub-ról
drugs <- 
  readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2023/2023-03-14/drugs.csv')

# Adatok áttekintése
dim(drugs)

## [1] 1988   28

names(drugs)

##  [1] "category"                                   
##  [2] "medicine_name"                              
##  [3] "therapeutic_area"                           
##  [4] "common_name"                                
##  [5] "active_substance"                           
##  [6] "product_number"                             
##  [7] "patient_safety"                             
##  [8] "authorisation_status"                       
##  [9] "atc_code"                                   
## [10] "additional_monitoring"                      
## [11] "generic"                                    
## [12] "biosimilar"                                 
## [13] "conditional_approval"                       
## [14] "exceptional_circumstances"                  
## [15] "accelerated_assessment"                     
## [16] "orphan_medicine"                            
## [17] "marketing_authorisation_date"               
## [18] "date_of_refusal_of_marketing_authorisation" 
## [19] "marketing_authorisation_holder_company_name"
## [20] "pharmacotherapeutic_group"                  
## [21] "date_of_opinion"                            
## [22] "decision_date"                              
## [23] "revision_number"                            
## [24] "condition_indication"                       
## [25] "species"                                    
## [26] "first_published"                            
## [27] "revision_date"                              
## [28] "url"

glimpse(drugs)

## Rows: 1,988
## Columns: 28
## $ category                                    <chr> "human", "human", "human",…
## $ medicine_name                               <chr> "Adcetris", "Nityr", "Ebva…
## $ therapeutic_area                            <chr> "Lymphoma, Non-Hodgkin;  H…
## $ common_name                                 <chr> "brentuximab vedotin", "ni…
## $ active_substance                            <chr> "brentuximab vedotin", "ni…
## $ product_number                              <chr> "002455", "004582", "00457…
## $ patient_safety                              <lgl> FALSE, FALSE, FALSE, FALSE…
## $ authorisation_status                        <chr> "authorised", "authorised"…
## $ atc_code                                    <chr> "L01XC12", "A16AX04", NA, …
## $ additional_monitoring                       <lgl> FALSE, FALSE, TRUE, TRUE, …
## $ generic                                     <lgl> FALSE, TRUE, FALSE, FALSE,…
## $ biosimilar                                  <lgl> FALSE, FALSE, FALSE, FALSE…
## $ conditional_approval                        <lgl> FALSE, FALSE, FALSE, FALSE…
## $ exceptional_circumstances                   <lgl> FALSE, FALSE, TRUE, FALSE,…
## $ accelerated_assessment                      <lgl> FALSE, FALSE, FALSE, FALSE…
## $ orphan_medicine                             <lgl> TRUE, FALSE, TRUE, FALSE, …
## $ marketing_authorisation_date                <date> 2012-10-25, 2018-07-26, 2…
## $ date_of_refusal_of_marketing_authorisation  <date> NA, NA, NA, NA, NA, NA, N…
## $ marketing_authorisation_holder_company_name <chr> "Takeda Pharma A/S", "Cycl…
## $ pharmacotherapeutic_group                   <chr> "Antineoplastic agents", "…
## $ date_of_opinion                             <date> 2012-07-19, 2018-05-31, 2…
## $ decision_date                               <date> 2022-11-17, 2023-03-10, 2…
## $ revision_number                             <dbl> 34, 4, 2, 3, 30, 24, 4, 18…
## $ condition_indication                        <chr> "Hodgkin lymphomaAdcetris …
## $ species                                     <chr> NA, NA, NA, NA, NA, NA, NA…
## $ first_published                             <dttm> 2018-07-25 13:58:00, 2018…
## $ revision_date                               <dttm> 2023-03-13 11:52:00, 2023…
## $ url                                         <chr> "https://www.ema.europa.eu…

skimr::skim(drugs)

Data summary
Name	drugs
Number of rows	1988
Number of columns	28
_______________________
Column type frequency:
character	13
Date	4
logical	8
numeric	1
POSIXct	2
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	n_unique
category	0	1.00	5	10	2
medicine_name	0	1.00	3	125	1976
therapeutic_area	285	0.86	4	400	669
common_name	4	1.00	4	220	1261
active_substance	1	1.00	4	823	1345
product_number	0	1.00	6	6	1932
authorisation_status	1	1.00	7	10	3
atc_code	28	0.99	3	18	1074
marketing_authorisation_holder_company_name	4	1.00	4	65	615
pharmacotherapeutic_group	34	0.98	7	174	365
condition_indication	12	0.99	18	7597	1886
species	1709	0.14	4	67	59
url	0	1.00	53	148	1988

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
marketing_authorisation_date	60	0.97	1995-10-20	2023-02-20	2013-06-09	1127
date_of_refusal_of_marketing_authorisation	1913	0.04	2004-09-07	2022-04-29	2013-04-25	67
date_of_opinion	779	0.61	1995-07-12	2022-12-15	2016-07-21	389
decision_date	45	0.98	1998-08-20	2023-03-10	2022-02-16	815

Variable type: logical

skim_variable	complete_rate	mean	count
patient_safety	1	0.01	FAL: 1977, TRU: 11
additional_monitoring	1	0.19	FAL: 1601, TRU: 387
generic	1	0.16	FAL: 1673, TRU: 315
biosimilar	1	0.05	FAL: 1896, TRU: 92
conditional_approval	1	0.02	FAL: 1940, TRU: 48
exceptional_circumstances	1	0.02	FAL: 1940, TRU: 48
accelerated_assessment	1	0.02	FAL: 1940, TRU: 48
orphan_medicine	1	0.08	FAL: 1826, TRU: 162

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
revision_number	96	0.95	13.53	11.65	0	4.75	11	19	89	▇▃▁▁▁

Variable type: POSIXct

skim_variable	n_missing	complete_rate	min	max	median	n_unique
first_published	0	1.00	1998-08-20 00:00:00	2023-03-09 18:50:00	2018-04-14 23:14:30	1760
revision_date	29	0.99	2000-07-17 02:00:00	2023-03-13 11:52:00	2022-05-11 11:58:00	1932

# Folytonos változók kialakítása
## Gyógyszer piacon töltött idejének (days_on_market változó) kiszámítása napokban (Az adatbázis feltöltési dátuma: 2023.03.13.)
drugs <- drugs %>%
  mutate(
    days_on_market = as.numeric(ymd("2023-03-13") - ymd(marketing_authorisation_date)))

## Ellenörzés - days_on_market
summary(drugs$days_on_market)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      21    1818    3564    3948    5655   10006      60

## Eloszlás - days_on_market
ggplot(drugs, aes(x = days_on_market)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  labs(x = "Napok száma a piacon", y = "Gyakoriság") +
  theme_minimal()

A „piacon töltött időt” (days_on_market változó) a forgalomba hozatali engedély dátuma és az adatbázis feltöltési dátuma (2023-03-13) közötti napok száma alapján képeztem, ezzel biztosítva, hogy a mutató ne függjön a futtatás aktuális dátumától.

Adatok feldolgozása

Minta

Az adatbázis összesen 1988 gyógyszert tartalmaz, amelyek közül 1706 (85.8%) humán és 282 (14.2%) állatgyógyászati készítmény. Az adatbázis 28 változót tartalmaz, beleértve a gyógyszer nevét, terápiás területét, hatóanyagát, ATC kódját, engedélyezési dátumát és státuszát, valamint olyan különböző jellemzőket a gyógyszerekkel kapcsolatban, mint pl. generikus, orphan, gyorsított eljárás stb.

A gyógyszerek engedélyezési dátuma 1995. október 20. és 2023. február 20. között oszlik meg (medián: 2013. június 9.). Az engedélyezési státusz tekintetében a gyógyszerek 79.1%-a (n = 1573) rendelkezik érvényes engedéllyel, 18.0% (n = 357) visszavont, és 2.9% (n = 57) elutasított státuszú.

A speciális státuszok tekintetében a gyógyszerek 15.8%-a (n = 315) generikus, 8.2%-a (n = 162) orphan gyógyszer és 2.4%-a (n = 48) részesült gyorsított eljárásban.

A gyógyszerek piacon töltött idejének eloszlása a teljes mintában széles tartományt fed le (Min: 21 nap; Mdn: 3564 nap; M: 3948 nap; Max: 10006 nap; NA: 60).

# Kategórikus változók gyakorisága és százalékos eloszlása
vars <- c(
  "category", "patient_safety", "authorisation_status", "generic",
  "conditional_approval", "exceptional_circumstances",
  "accelerated_assessment", "orphan_medicine")

cat_tables <- drugs %>%
  select(all_of(vars)) %>%
  lapply(function(x) {
    tab <- table(x, useNA = "ifany")
    pct <- prop.table(tab) * 100
    data.frame(
      value   = names(tab),
      n       = as.integer(tab),
      percent = round(as.numeric(pct), 1),
      row.names = NULL)})

cat_tables$category

##        value    n percent
## 1      human 1706    85.8
## 2 veterinary  282    14.2

cat_tables$patient_safety

##   value    n percent
## 1 FALSE 1977    99.4
## 2  TRUE   11     0.6

cat_tables$authorisation_status

##        value    n percent
## 1 authorised 1573    79.1
## 2    refused   57     2.9
## 3  withdrawn  357    18.0
## 4       <NA>    1     0.1

cat_tables$generic

##   value    n percent
## 1 FALSE 1673    84.2
## 2  TRUE  315    15.8

cat_tables$conditional_approval

##   value    n percent
## 1 FALSE 1940    97.6
## 2  TRUE   48     2.4

cat_tables$exceptional_circumstances

##   value    n percent
## 1 FALSE 1940    97.6
## 2  TRUE   48     2.4

cat_tables$accelerated_assessment

##   value    n percent
## 1 FALSE 1940    97.6
## 2  TRUE   48     2.4

cat_tables$orphan_medicine

##   value    n percent
## 1 FALSE 1826    91.9
## 2  TRUE  162     8.1

# Piacon töltött idő (days_on_market) leíró statisztikája
summary(drugs$days_on_market)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      21    1818    3564    3948    5655   10006      60

# Gyógyszerek engedélyezésének száma évente
drugs %>%
  filter(!is.na(marketing_authorisation_date)) %>%
  mutate(year = year(marketing_authorisation_date)) %>%
  count(year) %>%
  ggplot(aes(year, n)) +
  geom_line() +
  scale_x_continuous(breaks = seq(1995, 2022, by = 1), minor_breaks = NULL, limits = c(NA, 2022)) +
  labs(title = "Gyógyszer-engedélyezések száma évente", x = "Év", y = "n") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Top 10 gyártó
drugs %>%
  count(marketing_authorisation_holder_company_name, sort = TRUE) %>%
  head(10)

## # A tibble: 10 × 2
##    marketing_authorisation_holder_company_name     n
##    <chr>                                       <int>
##  1 Accord Healthcare S.L.U.                       58
##  2 Novartis Europharm Limited                     58
##  3 Pfizer Europe MA EEIG                          43
##  4 Zoetis Belgium SA                              40
##  5 AstraZeneca AB                                 36
##  6 Boehringer Ingelheim Vetmedica GmbH            36
##  7 Merck Sharp & Dohme B.V.                       32
##  8 Teva B.V.                                      32
##  9 Intervet International BV                      31
## 10 Eli Lilly Nederland B.V.                       30

Gyógyszer-engedélyezési trendek

Az 1. ábra a gyógyszer-engedélyezések számának alakulását mutatja az 1995-2022 közötti időszakban. A trendvonal jelentős fluktuációt mutat, kiugró csúcsokkal 2001 (n ≈ 90), 2017 (n ≈ 110) és 2021 körül (n ≈ 120). A COVID-19 pandémia időszakában (2020-2021) a gyógyszer-engedélyezések száma jelentősen megugrott, ami feltételezhetően a vakcinák és vírusellenes szerek gyorsított engedélyezésének tudható be elsősorban.

A forgalmazási engedély jogosultjai közül a legtöbb gyógyszert az Accord Healthcare S.L.U. és a Novartis Europharm Limited jegyzi (egyaránt n = 58), amelyeket a Pfizer Europe MA EEIG (n = 43) és a Zoetis Belgium SA (n = 40) követ.

ATC-rendszer

A humán készítmények ATC-főcsoport szerinti bontása lehetővé teszi annak áttekintését, hogy az engedélyezések elsősorban mely terápiás területekhez kapcsolódtak, illetve hogyan alakultak ezen időszakban. Az ATC-főcsoportok szerinti évekre vontott ábra a főcsoportok közötti nagyságrendi különbségeket és az időbeli dinamikát együttesen teszi láthatóvá, így a későbbi, célzottabb összehasonlítások megalapozására szolgálhat.

# Human gyógyszerek ATC kód csoportonként
drugs %>%
  filter(category == "human", !is.na(marketing_authorisation_date), !is.na(atc_code)) %>%
  mutate(
    year = year(marketing_authorisation_date),
    atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  count(year, atc_group) %>%
  ggplot(aes(year, n)) +
  geom_line() +
  facet_wrap(~atc_group, scales = "free_x") +
  scale_x_continuous(breaks = seq(1995, 2022, by = 5), limits = c(1995, 2022), minor_breaks = NULL) +
  labs(
    title = "Humán gyógyszer-engedélyezések száma évente ATC-rendszer szerint",
    x = "Év",
    y = "n",
    caption = "Note. ATC besorolás: A - Tápcsatorna és anyagcsere; B - Vér és vérképzőszervek; C - Cardiovascularis rendszer; D - Bőrgyógyászati készítmények; G - Urogenitalis rendszer és nemi hormonok;\nH - Szisztémás hormonkészítmények; J - Szisztémás fertőzésellenes szerek; L - Daganatellenes és immunmoduláns szerek; M - Váz- és izomrendszer; N - Központi idegrendszer;\nP - Parazitaellenes készítmények; R - Légzőrendszer; S - Érzékszervek; V - Egyéb") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text(size = 10),
    plot.caption = element_text(hjust = 0, size = 9))

# Heatmap - ATC csoportok évenete
drugs %>%
  filter(category == "human", !is.na(marketing_authorisation_date), !is.na(atc_code)) %>%
  mutate(
    year = year(marketing_authorisation_date),
    atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  filter(year <= 2022) %>%
  count(year, atc_group) %>%
  group_by(year) %>%
  mutate(percent = n / sum(n) * 100) %>%
  ungroup() %>%
  ggplot(aes(x = atc_group, y = factor(year), fill = percent)) +
  geom_tile(color = "white") +
  geom_text(aes(label = ifelse(n > 0, n, "")), size = 2.5) +
  scale_fill_gradient(low = "white", high = "steelblue", na.value = "white") +
  labs(
    title = "Humán gyógyszer-engedélyezésének száma ATC-rendszer szerint",
    x = "ATC csoport",
    y = "Év",
    fill = "Arány (%)",
    caption = "Note. ATC besorolás: A - Tápcsatorna és anyagcsere; B - Vér és vérképzőszervek; C - Cardiovascularis rendszer; D - Bőrgyógyászati készítmények; G - Urogenitalis rendszer és nemi hormonok; \nH - Szisztémás hormonkészítmények; J - Szisztémás fertőzésellenes szerek; L - Daganatellenes és immunmoduláns szerek; M - Váz- és izomrendszer; N - Központi idegrendszer; \nP - Parazitaellenes készítmények; R - Légzőrendszer; S - Érzékszervek; V - Egyéb") +
  theme_minimal() +
  theme(
    plot.caption = element_text(hjust = 0, size = 9),
    panel.grid = element_blank())

# ATC csoportok gyakorisága a teljes időszakban
drugs %>%
  filter(category == "human", !is.na(atc_code)) %>%
  mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  count(atc_group, sort = TRUE) %>%
  mutate(percent = round(n / sum(n) * 100, 1))

## # A tibble: 14 × 3
##    atc_group     n percent
##    <chr>     <int>   <dbl>
##  1 L           463    27.5
##  2 J           253    15.1
##  3 A           201    12  
##  4 N           186    11.1
##  5 B           135     8  
##  6 C            95     5.7
##  7 V            72     4.3
##  8 R            69     4.1
##  9 G            59     3.5
## 10 M            54     3.2
## 11 H            40     2.4
## 12 S            34     2  
## 13 D            18     1.1
## 14 P             2     0.1

Az ATC-rendszer szerinti megoszlás

Az ATC-rendszer szerinti elemzés a humán gyógyszerekre korlátozódott. A 14 fő ATC csoportból a leggyakoribb az L csoport (daganatellenes és immunmoduláns szerek; n = 463, 27.5%), amelyet a J csoport (szisztémás fertőzésellenes szerek; n = 253, 15.1%), az A csoport (tápcsatorna és anyagcsere; n = 201, 12.0%) és az N csoport (központi idegrendszer; n = 186, 11.1%) követ. Ez a négy csoport együttesen a humán gyógyszerek 65.7%-át teszi ki.

A heatmap-vizualizáció (2. ábra) az ATC csoportok éves megoszlását szemlélteti. Az L csoport dominanciája különösen az elmúlt évtizedben vált hangsúlyossá, ami feltehetően a precíziós onkológia és az immunterápiák fejlődésével magyarázható.

A további elemzéseket a négy leggyakoribb ATC csoportra (L, J, A, N) szűkítettem, amelyek együttesen 1103 gyógyszert foglalnak magukban.

Statisztikai elemzések - TOP4 ATC

# Top 4 ATC csoport kiszűrése
top4_atc <- drugs %>%
  filter(category == "human", !is.na(atc_code)) %>%
  mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  filter(atc_group %in% c("L", "J", "A", "N"))

# Leíró statisztika - revision_number
top4_atc %>%
  group_by(atc_group) %>%
  summarise(
    n = n(),
    hianyzo = sum(is.na(revision_number)),
    atlag = mean(revision_number, na.rm = TRUE),
    sd = sd(revision_number, na.rm = TRUE),
    min = min(revision_number, na.rm = TRUE),
    q25 = quantile(revision_number, 0.25, na.rm = TRUE),
    median = median(revision_number, na.rm = TRUE),
    q75 = quantile(revision_number, 0.75, na.rm = TRUE),
    max = max(revision_number, na.rm = TRUE))

## # A tibble: 4 × 10
##   atc_group     n hianyzo atlag    sd   min   q25 median   q75   max
##   <chr>     <int>   <int> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 A           201       7  14.1  9.51     0     6     14  21      37
## 2 J           253       6  17.7 14.2      0     7     14  26.5    68
## 3 L           463      31  14.5 12.9      0     5     11  20      89
## 4 N           186       8  15.0 11.6      0     7     13  21      58

# Leíró statisztika - days_on_market
top4_atc %>%
  group_by(atc_group) %>%
  summarise(
    n = n(),
    hianyzo = sum(is.na(days_on_market)),
    atlag = mean(days_on_market, na.rm = TRUE),
    sd = sd(days_on_market, na.rm = TRUE),
    min = min(days_on_market, na.rm = TRUE),
    q25 = quantile(days_on_market, 0.25, na.rm = TRUE),
    median = median(days_on_market, na.rm = TRUE),
    q75 = quantile(days_on_market, 0.75, na.rm = TRUE),
    max = max(days_on_market, na.rm = TRUE))

## # A tibble: 4 × 10
##   atc_group     n hianyzo atlag    sd   min   q25 median   q75   max
##   <chr>     <int>   <int> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 A           201       4 4182. 2524.    94 2163   3979  5908   9813
## 2 J           253       2 4300. 2774.    98 2008.  4042  6208.  9805
## 3 L           463      20 3130. 2455.    26 1138   2482  4676   9968
## 4 N           186       8 3987. 2420.   126 2267   3808. 5258.  9772

# ATC csoport nevek hozzáadása
top4_atc <- top4_atc %>%
  mutate(atc_name = case_when(
    atc_group == "L" ~ "L - Daganatellenes és\nimmunomoduláns szerek",
    atc_group == "J" ~ "J - Szisztémás\nfertőzésellenes szerek",
    atc_group == "A" ~ "A - Tápcsatorna\nés anyagcsere",
    atc_group == "N" ~ "N - Központi\nidegrendszer"
  ))

# Violin plot - revision_number
ggplot(top4_atc, aes(x = atc_name, y = revision_number, fill = atc_name)) +
  geom_violin(alpha = 0.5) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
  labs(
    title = "Revíziók száma ATC csoportonként",
    x = "ATC csoport",
    y = "Revíziók száma"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Violin plot - days_on_market
ggplot(top4_atc, aes(x = atc_name, y = days_on_market, fill = atc_name)) +
  geom_violin(alpha = 0.5) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
  labs(
    title = "Piacon eltöltött napok száma ATC csoportonként",
    x = "ATC csoport",
    y = "Napok száma"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

Leíró statisztika

A revíziók számának leíró statisztikáját az 1. táblázat foglalja össze. A revíziók átlagos száma a négy ATC csoportban 14.1 és 17.7 között mozgott. A legmagasabb átlagot a J csoport mutatta (M = 17.7, SD = 14.2), míg a legalacsonyabbat az A csoport (M = 14.1, SD = 9.5).

A piacon töltött idő tekintetében az L csoport gyógyszerei átlagosan rövidebb ideje vannak a piacon (M = 3130 nap, SD = 2455), mint az A (M = 4182 nap, SD = 2524) és J (M = 4300 nap, SD = 2774) csoport készítményei.

A violin plotok (3. és 4. ábra) vizuálisan is szemléltetik a csoportok közötti különbségeket.

# Függő változó: revision_number
## Shapiro-Wilk teszt - revision_number
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    W = round(shapiro.test(revision_number)$statistic, 3),
    p_value = format_p(shapiro.test(revision_number)$p.value))

## # A tibble: 4 × 3
##   atc_name                                           W p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.959 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.922 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.87  < 0.001
## 4 "N - Központi\nidegrendszer"                   0.919 < 0.001

## Lilliefors teszt - revision_number
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    D = round(lillie.test(revision_number)$statistic, 3),
    p_value = format_p(lillie.test(revision_number)$p.value))

## # A tibble: 4 × 3
##   atc_name                                           D p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.108 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.138 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.132 < 0.001
## 4 "N - Központi\nidegrendszer"                   0.123 < 0.001

# Függő változó: days_on_market
## Shapiro-Wilk teszt - days_on_market
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    W = round(shapiro.test(days_on_market)$statistic, 3),
    p_value = format_p(shapiro.test(days_on_market)$p.value))

## # A tibble: 4 × 3
##   atc_name                                           W p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.965 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.945 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.911 < 0.001
## 4 "N - Központi\nidegrendszer"                   0.958 < 0.001

## Lilliefors teszt - days_on_market
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    D = round(lillie.test(days_on_market)$statistic, 3),
    p_value = format_p(lillie.test(days_on_market)$p.value))

## # A tibble: 4 × 3
##   atc_name                                           D p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.069 0.025  
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.094 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.116 < 0.001
## 4 "N - Központi\nidegrendszer"                   0.062 0.091

# Ferdeség és csúcsosság
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    skewness_revision = round(skewness(revision_number, na.rm = TRUE), 3),
    kurtosis_revision = round(kurtosis(revision_number, na.rm = TRUE), 3),
    skewness_days = round(skewness(days_on_market, na.rm = TRUE), 3),
    kurtosis_days = round(kurtosis(days_on_market, na.rm = TRUE), 3))

## # A tibble: 4 × 5
##   atc_name       skewness_revision kurtosis_revision skewness_days kurtosis_days
##   <chr>                      <dbl>             <dbl>         <dbl>         <dbl>
## 1 "A - Tápcsato…             0.357              2.19         0.286          2.12
## 2 "J - Szisztém…             0.93               3.34         0.378          2.02
## 3 "L - Daganate…             1.55               6.48         0.912          3.05
## 4 "N - Központi…             1.07               4.00         0.528          2.73

# Vizuális módszerek
## Hisztogram - revision_number
ggplot(top4_atc, aes(x = revision_number, fill = atc_name)) +
  geom_histogram(bins = 30, color = "white") +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Revíziók számának eloszlása ATC csoportonként",
    x = "Revíziók száma",
    y = "Gyakoriság") +
  theme_minimal() +
  theme(legend.position = "none")

## Q-Q plot - revision_number
ggplot(top4_atc, aes(sample = revision_number, color = atc_name)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Q-Q plot: Revíziók száma ATC csoportonként",
    x = "Elméleti kvantilisek",
    y = "Minta kvantilisek") +
  theme_minimal() +
  theme(legend.position = "none")

## Hisztogram - days_on_market
ggplot(top4_atc, aes(x = days_on_market, fill = atc_name)) +
  geom_histogram(bins = 30, color = "white") +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Piacon eltöltött napok eloszlása ATC csoportonként",
    x = "Napok száma",
    y = "Gyakoriság") +
  theme_minimal() +
  theme(legend.position = "none")

## Q-Q plot - days_on_market
ggplot(top4_atc, aes(sample = days_on_market, color = atc_name)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Q-Q plot: Piacon eltöltött napok ATC csoportonként",
    x = "Elméleti kvantilisek",
    y = "Minta kvantilisek") +
  theme_minimal() +
  theme(legend.position = "none")

A normalitás vizsgálata

A paraméteres statisztikai eljárások alkalmazhatóságának előfeltételeként a függő változók normalitását Shapiro-Wilk és Lilliefors (Kolmogorov-Smirnov korrekció) tesztekkel vizsgáltam.

A revíziók számának eloszlása a Shapiro-Wilk teszt alapján egyik ATC csoportban sem követte a normális eloszlást (p < .001). A Lilliefors teszt eredményei megerősítették ezt a mintázatot (p < .001). A ferdeségi mutatók pozitív értékeket vettek fel (skewness = 0.36–1.55), ami jobbra ferde eloszlásokra utal. A csúcsossági értékek (kurtosis = 2.19–6.48) az L csoportban kifejezetten leptokurtikus eloszlást jeleztek.

A piacon töltött idő esetében hasonló eredmények lárhatók. A Shapiro-Wilk teszt valamennyi csoportban szignifikáns eltérést mutatott a normalitástól (p < .001). A Lilliefors teszt az N csoport kivételével (p = .091) szintén szignifikáns eredményeket adott.

A Q-Q plotok (5-8. ábra) vizuálisan is alátámasztják a normalitástól való eltérést, különösen az eloszlások szélein megfigyelhető szisztematikus eltérésekkel.

Tekintettel a normalitás sérülésére, a csoportok összehasonlításához nem-paraméteres statisztikai eljárásokat alkalmaztam.

# Függő változó: revison_number
## Kruskal-Wallis teszt - revision_number
top4_atc %>%
  kruskal_test(revision_number ~ atc_group) %>%
  mutate(p = format_p(p))

## # A tibble: 1 × 6
##   .y.                 n statistic    df p     method        
##   <chr>           <int>     <dbl> <int> <chr> <chr>         
## 1 revision_number  1103      9.39     3 0.025 Kruskal-Wallis

## Kruskal-Wallis hatásméret (eta²) - revision_number
top4_atc %>%
  kruskal_effsize(revision_number ~ atc_group, ci = TRUE)

## # A tibble: 1 × 7
##   .y.                 n effsize conf.low conf.high method  magnitude
## * <chr>           <int>   <dbl>    <dbl>     <dbl> <chr>   <ord>    
## 1 revision_number  1103 0.00581 -0.00073      0.02 eta2[H] small

## Dunn-teszt Holm korrekcióval - revision_number
dunn_revision <- dunn_test(revision_number ~ atc_group, data = top4_atc, p.adjust.method = "holm")
dunn_revision %>%
  mutate(p.adj = format_p(p.adj))

## # A tibble: 6 × 9
##   .y.             group1 group2    n1    n2 statistic       p p.adj p.adj.signif
##   <chr>           <chr>  <chr>  <int> <int>     <dbl>   <dbl> <chr> <chr>       
## 1 revision_number A      J        194   247    1.60   0.110   0.552 ns          
## 2 revision_number A      L        194   432   -1.05   0.292   0.833 ns          
## 3 revision_number A      N        194   178    0.0537 0.957   0.957 ns          
## 4 revision_number J      L        247   432   -3.06   0.00219 0.013 *           
## 5 revision_number J      N        247   178   -1.50   0.133   0.552 ns          
## 6 revision_number L      N        432   178    1.09   0.278   0.833 ns

## Cliff's delta függvény 95% CI-vel
calc_cliffs_delta <- function(data, var, group1, group2) {
  x <- data %>% filter(atc_group == group1) %>% pull({{var}})
  y <- data %>% filter(atc_group == group2) %>% pull({{var}})
  result <- effectsize::cliffs_delta(x, y, ci = 0.95)
  tibble(
    group1 = group1,
    group2 = group2,
    cliffs_delta = round(result$r, 3),
    ci_low = round(result$CI_low, 3),
    ci_high = round(result$CI_high, 3)
  )
}

## Cliff's delta páronként - revision_number
pairs <- combn(unique(top4_atc$atc_group), 2, simplify = FALSE)
map_df(pairs, ~calc_cliffs_delta(top4_atc, revision_number, .x[1], .x[2]))

## # A tibble: 6 × 5
##   group1 group2 cliffs_delta ci_low ci_high
##   <chr>  <chr>         <dbl>  <dbl>   <dbl>
## 1 L      A            -0.055 -0.152   0.043
## 2 L      J            -0.137 -0.224  -0.048
## 3 L      N            -0.059 -0.158   0.042
## 4 A      J            -0.094 -0.201   0.014
## 5 A      N            -0.001 -0.118   0.116
## 6 J      N             0.088 -0.023   0.197

# Függő változó: days_on_market
## Kruskal-Wallis teszt - days_on_market
top4_atc %>%
  kruskal_test(days_on_market ~ atc_group) %>%
  mutate(p = format_p(p))

## # A tibble: 1 × 6
##   .y.                n statistic    df p       method        
##   <chr>          <int>     <dbl> <int> <chr>   <chr>         
## 1 days_on_market  1103      48.2     3 < 0.001 Kruskal-Wallis

## Kruskal-Wallis hatásméret (eta²) - days_on_market
top4_atc %>%
  kruskal_effsize(days_on_market ~ atc_group, ci = TRUE)

## # A tibble: 1 × 7
##   .y.                n effsize conf.low conf.high method  magnitude
## * <chr>          <int>   <dbl>    <dbl>     <dbl> <chr>   <ord>    
## 1 days_on_market  1103  0.0412     0.02      0.07 eta2[H] small

## Dunn-teszt Holm korrekcióval - days_on_market
dunn_days <- dunn_test(days_on_market ~ atc_group, data = top4_atc, p.adjust.method = "holm")
dunn_days %>%
  mutate(p.adj = format_p(p.adj))

## # A tibble: 6 × 9
##   .y.            group1 group2    n1    n2 statistic        p p.adj p.adj.signif
##   <chr>          <chr>  <chr>  <int> <int>     <dbl>    <dbl> <chr> <chr>       
## 1 days_on_market A      J        197   251     0.105  9.17e-1 1     ns          
## 2 days_on_market A      L        197   443    -5.15   2.58e-7 < 0.… ****        
## 3 days_on_market A      N        197   178    -0.551  5.82e-1 1     ns          
## 4 days_on_market J      L        251   443    -5.71   1.13e-8 < 0.… ****        
## 5 days_on_market J      N        251   178    -0.683  4.94e-1 1     ns          
## 6 days_on_market L      N        443   178     4.33   1.50e-5 < 0.… ****

## Cliff's delta páronként - days_on_market
map_df(pairs, ~calc_cliffs_delta(top4_atc, days_on_market, .x[1], .x[2]))

## # A tibble: 6 × 5
##   group1 group2 cliffs_delta ci_low ci_high
##   <chr>  <chr>         <dbl>  <dbl>   <dbl>
## 1 L      A            -0.256 -0.344  -0.163
## 2 L      J            -0.254 -0.336  -0.169
## 3 L      N            -0.23  -0.323  -0.133
## 4 A      J            -0.013 -0.12    0.095
## 5 A      N             0.04  -0.077   0.156
## 6 J      N             0.047 -0.064   0.157

Csoportok összehasonlítása

Revíziók száma

A négy ATC csoport közötti különbséget Kruskal-Wallis teszttel vizsgáltam. Az eredmények szignifikáns különbséget mutattak a csoportok között (H(3) = 9.39, p = .025). A hatásméret ugyanakkor kicsi volt (η² = 0.006, 95% CI [-0.001, 0.020]). A 95% CI alsó határa 0 közelébe esik (a becslési eljárás miatt numerikusan enyhén negatív), ami a hatás elhanyagolható mértékével konzisztens. A Dunn-féle post-hoc teszt Holm-korrekcióval egyetlen szignifikáns páronkénti különbséget tárt fel. A J csoport (szisztémás fertőzésellenes szerek) szignifikánsan magasabb revíziószámmal rendelkezett, mint az L csoport (daganatellenes szerek), (padj = .013). A Cliff’s delta hatásméret mutató eredménye alapján kis hatásméret írható le a két csoport között (d = -0.137, 95% CI [-0.224, -0.048]).

Piacon töltött idő

A piacon töltött idő tekintetében a Kruskal-Wallis teszt erősen szignifikáns különbséget mutatott a csoportok között (H(3) = 48.2, p < .001). A hatásméret kis-közepes mértékű volt (η² = 0.041, 95% CI [0.020, 0.070]). A Dunn-féle post-hoc teszt három szignifikáns páronkénti különbséget tárt fel (padj < .001). Az L csoport gyógyszerei szignifikánsan rövidebb ideje vannak a piacon, mint az A csoport (d = -0.256), a J csoport (d = -0.254) és az N csoport (d = -0.230) készítményei. Ez a mintázat összhangban van az onkológiai területen tapasztalható intenzívebb kutatás-fejlesztési aktivitással és a gyorsabb gyógyszercserélődéssel.

Hálózatelemzés

Adatok előkészítése

# A therapeutic_area változó tördelése
drugs_clean <- drugs %>%
  separate_rows(therapeutic_area, sep = "\\s*;\\s*") %>% 
  mutate(therapeutic_area = trimws(therapeutic_area))

# N ATC csoport kiszűrése (központi idegrendszer)
drugs_N <- drugs_clean %>%
  filter(!is.na(atc_code)) %>%
  mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  filter(atc_group == "N") %>%
  filter(!is.na(therapeutic_area))

# Összefoglaló statisztika
drugs_N %>%
  summarise(
    n_drugs = n_distinct(medicine_name),
    n_areas = n_distinct(therapeutic_area),
    n_substances = n_distinct(active_substance),
    n_companies = n_distinct(marketing_authorisation_holder_company_name))

## # A tibble: 1 × 4
##   n_drugs n_areas n_substances n_companies
##     <int>   <int>        <int>       <int>
## 1     185      70          101         111

A hálózatelemzés a központi idegrendszeri (N) ATC csoportra fókuszált, amely a pszichiátriai és neurológiai gyógyszereket foglalja magában. A vizsgálat célja a terápiás területek és hatóanyagok közötti kapcsolatok feltérképezése volt.

Az N csoportban összesen 185 egyedi gyógyszer, 70 terápiás terület, 101 hatóanyag és 111 forgalmazó cég volt azonosítható. A therapeutic_area változó pontosvesszővel elválasztott többszörös értékeket tartalmazott, amelyeket külön sorokra bontottam a hálózatépítéséhez.

Szűrés és hálózat építése

# Csak a gyakoribb terápiás területek és hatóanyagok (min. 3 előfordulás)
top_areas_N <- drugs_N %>%
  count(therapeutic_area) %>%
  filter(n >= 3) %>%
  pull(therapeutic_area)

top_substances_N <- drugs_N %>%
  count(active_substance) %>%
  filter(n >= 3) %>%
  pull(active_substance)

# Szűrt adatok
drugs_N_filtered <- drugs_N %>%
  filter(therapeutic_area %in% top_areas_N | active_substance %in% top_substances_N)

# Élsúlyok hozzáadása (hány gyógyszer kapcsolja össze őket)
edges_weighted <- drugs_N_filtered %>%
  count(therapeutic_area, active_substance, name = "weight") %>%
  rename(from = therapeutic_area, to = active_substance)

# ATC alcsoport kinyerése (N01, N02, stb.)
atc_subgroups <- drugs_N_filtered %>%
  select(therapeutic_area, active_substance, atc_code) %>%
  mutate(atc_subgroup = str_extract(atc_code, "^N\\d{2}")) %>%
  distinct(active_substance, atc_subgroup) %>%
  group_by(active_substance) %>%
  slice(1) %>%
  ungroup()

# Csúcsok létrehozása típussal és ATC alcsoporttal
nodes_areas_f <- drugs_N_filtered %>%
  distinct(therapeutic_area) %>%
  rename(name = therapeutic_area) %>%
  mutate(type = "Terápiás terület", atc_subgroup = NA_character_)

nodes_substances_f <- drugs_N_filtered %>%
  distinct(active_substance) %>%
  rename(name = active_substance) %>%
  mutate(type = "Hatóanyag") %>%
  left_join(atc_subgroups, by = c("name" = "active_substance"))

nodes_f <- bind_rows(nodes_areas_f, nodes_substances_f)

# Hálózat létrehozása súlyozott élekkel
network_f <- tbl_graph(nodes = nodes_f, edges = edges_weighted, directed = FALSE)

# Fokszám és klaszter hozzáadása
network_f <- network_f %>%
  activate(nodes) %>%
  mutate(
    degree = centrality_degree(),
    cluster = as.factor(group_louvain()))

Hálózat vizualizáció

# Hálózat vizualizációja - Interaktív ábra (jobb átláthatóság)
layout_coords <- create_layout(network_f, layout = "fr")
layout_coords$x <- layout_coords$x * 2
layout_coords$y <- layout_coords$y * 2

p <- ggraph(layout_coords) +
  geom_edge_link(aes(width = weight), alpha = 0.3, color = "gray50") +
  geom_point_interactive(aes(
    x = x, y = y,
    size = degree,
    color = ifelse(type == "Terápiás terület", "Terápiás terület", atc_subgroup),
    alpha = degree,
    tooltip = paste0(name, "\nKapcsolatok: ", degree),
    data_id = name
  )) +
  scale_size_continuous(range = c(2, 12)) +
  scale_alpha_continuous(range = c(0.5, 1)) +
  scale_edge_width_continuous(range = c(0.3, 3)) +
  scale_color_manual(
    values = c(
      "Terápiás terület" = "steelblue",
      "N01" = "#E41A1C",
      "N02" = "#FF7F00",
      "N03" = "#4DAF4A",
      "N04" = "#984EA3",
      "N05" = "#F781BF",
      "N06" = "#A65628",
      "N07" = "#999999"
    ),
    na.value = "coral"
  ) +
  labs(
    title = "Terápiás területek és hatóanyagok hálózata (Interaktív)",
    subtitle = "N (központi idegrendszer) ATC csoport - Kérem, vigye az egeret a pontokra a részletekért",
    color = "Típus / ATC alcsoport",
    size = "Kapcsolatok száma",
    edge_width = "Közös gyógyszerek"
  ) +
  theme_graph() +
  theme(legend.position = "right") +
  guides(alpha = "none")

girafe(
  ggobj = p,
  width_svg = 14,
  height_svg = 10,
  options = list(
    opts_zoom(max = 5),
    opts_hover(css = "fill:red;stroke:black;stroke-width:2;"),
    opts_tooltip(css = "background-color:white;padding:8px;border-radius:5px;border:1px solid gray;font-size:12px;")))

# Alapvető hálózati statisztikák
network_stats <- tibble(
  Mutató = c(
    "Csúcsok száma (nodes)",
    "Élek száma (edges)", 
    "Sűrűség (density)",
    "Átlagos fokszám (mean degree)",
    "Átmérő (diameter)",
    "Átlagos úthossz (mean distance)",
    "Tranzitivitás (clustering coefficient)"
  ),
  Érték = c(
    vcount(network_f),
    ecount(network_f),
    round(edge_density(network_f), 3),
    round(mean(degree(network_f)), 2),
    diameter(network_f),
    round(mean_distance(network_f), 2),
    round(transitivity(network_f), 3)))
network_stats

## # A tibble: 7 × 2
##   Mutató                                   Érték
##   <chr>                                    <dbl>
## 1 Csúcsok száma (nodes)                  110    
## 2 Élek száma (edges)                     108    
## 3 Sűrűség (density)                        0.018
## 4 Átlagos fokszám (mean degree)            1.96 
## 5 Átmérő (diameter)                       32    
## 6 Átlagos úthossz (mean distance)          6.61 
## 7 Tranzitivitás (clustering coefficient)   0

# Centralitás mutatók hozzáadása a hálózathoz
network_f <- network_f %>%
  activate(nodes) %>%
  mutate(
    degree = centrality_degree(),
    betweenness = round(centrality_betweenness(), 2),
    closeness = round(centrality_closeness(), 4),
    eigenvector = round(centrality_eigen(), 3))

# Top 20 csúcs különböző centralitás mutatók alapján
network_f %>%
  activate(nodes) %>%
  as_tibble() %>%
  arrange(desc(degree)) %>%
  select(name, type, degree, betweenness, closeness, eigenvector) %>%
  head(20)

## # A tibble: 20 × 6
##    name                           type  degree betweenness closeness eigenvector
##    <chr>                          <chr>  <dbl>       <dbl>     <dbl>       <dbl>
##  1 Parkinson Disease              Terá…     14       167.     0.0345       1    
##  2 amifampridine phosphate        Ható…     13       150      0.0312       0    
##  3 Schizophrenia                  Terá…     11        69.5    0.0625       0    
##  4 Epilepsy                       Terá…     10       176.     0.02         0    
##  5 Migraine Disorders             Terá…      6        15      0.167        0    
##  6 Alzheimer Disease              Terá…      5        54.3    0.0204       0.234
##  7 Opioid-Related Disorders       Terá…      5        10      0.2          0    
##  8 Pain                           Terá…      5        51.5    0.0172       0    
##  9 duloxetine                     Ható…      5        93      0.0145       0    
## 10 Bipolar Disorder               Terá…      4        13.5    0.0333       0    
## 11 Epilepsies, Partial            Terá…      4        41.5    0.0125       0    
## 12 sodium oxybate                 Ható…      4        12      0.125        0    
## 13 Sleep Initiation and Maintena… Terá…      3         5      0.2          0    
## 14 Depressive Disorder, Major     Terá…      3        41      0.0116       0    
## 15 Cancer                         Terá…      3        78.5    0.0263       0    
## 16 Neuralgia                      Terá…      3        63      0.0167       0    
## 17 Narcolepsy                     Terá…      3         9      0.111        0    
## 18 Pain, Postoperative            Terá…      3         3      0.333        0    
## 19 rivastigmine                   Ható…      3        37      0.0263       0.365
## 20 pregabalin                     Ható…      3       118      0.0189       0

Hálózatemezés mutatók

A kétoldalú (bipartite) hálózat 110 csúcspontot (terápiás területek és hatóanyagok) és 108 élt (kapcsolatot) tartalmaz. A hálózat ritkának tekinthető (sűrűség = 0.018), ami azt jelenti, hogy az összes lehetséges kapcsolat mindössze 1.8%-a realizálódik. Az átlagos fokszám 1.96, ami arra utal, hogy egy csúcspont átlagosan közel két másik csúcsponttal áll kapcsolatban.

A hálózat átmérője 32, az átlagos úthossz pedig 6.61, ami viszonylag szétszórt, kevéssé integrált struktúrára utal. A tranzitivitás (klaszterezési együttható) értéke 0, ami a bipartite hálózatok sajátossága, ahol azonos típusú csúcspontok nem kapcsolódnak közvetlenül egymáshoz.

Centralitás elemzés

A fokszám-centralitás alapján a legkiemelkedőbb terápiás területek a Parkinson-kór (fokszám = 14), a Szkizofrénia (fokszám = 11) és az Epilepszia (fokszám = 10). Ezek a betegségterületek számos különböző hatóanyaggal állnak kapcsolatban, ami a terápiás lehetőségek sokféleségét tükrözi.

A közöttiség-centralitás (betweenness) tekintetében az Epilepszia (176.0) és a Parkinson-kór (167.0) kiemelkedik, ami arra utal, hogy ezek a terápiás területek híd szerepet töltenek be a hálózatban, összekapcsolva más, egyébként elszigetelt klasztereket.

A hatóanyagok közül az amifampridine-foszfát rendelkezik a legmagasabb fokszámmal (13), amely a Lambert-Eaton myastheniás szindróma és más neuromuszkuláris betegségek kezelésében használatos. A duloxetin (fokszám = 5) és a pregabalin (fokszám = 3) szintén kiemelkedő pozíciót foglal el, ami széles terápiás alkalmazhatóságukat tükrözi (fájdalom, depresszió, szorongás, neuropátia).

Regressziós modellek

# Regressziós modell előkészítése
model_df <- top4_atc %>%
  transmute(
    revision_number,
    days_on_market,
    atc_group,
    generic,
    orphan_medicine,
    accelerated_assessment
  ) %>%
  filter(
    !is.na(revision_number),
    !is.na(days_on_market),
    !is.na(atc_group)
  ) %>%
  mutate(
    atc_group = factor(atc_group),
    generic = factor(generic),
    orphan_medicine = factor(orphan_medicine),
    accelerated_assessment = factor(accelerated_assessment),
    days_on_market_z = as.numeric(scale(days_on_market)))

skimr::skim(model_df)

Data summary
Name	model_df
Number of rows	1026
Number of columns	7
_______________________
Column type frequency:
factor	4
numeric	3
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
atc_group	1	FALSE	4	L: 418, J: 245, A: 190, N: 173
generic	1	FALSE	2	FAL: 837, TRU: 189
orphan_medicine	1	FALSE	2	FAL: 931, TRU: 95
accelerated_assessment	1	FALSE	2	FAL: 996, TRU: 30

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
revision_number	1	15.62	12.45	0.00	6.00	13.00	22.00	89.00	▇▃▁▁▁
days_on_market	1	3886.63	2542.55	124.00	1816.50	3462.50	5488.25	9968.00	▇▇▆▂▃
days_on_market_z	1	0.00	1.00	-1.48	-0.81	-0.17	0.63	2.39	▇▇▆▂▃

A revíziók számát befolyásoló tényezők azonosítására többváltozós lineáris regressziós modelleket építettem. A függő változó a revíziók száma volt, a prediktorok pedig az ATC csoport, a piacon töltött idő (z-standardizált), a generikus státusz, az orphan gyógyszer státusz és a gyorsított értékelés.

Az elemzés 1026 megfigyelésen alapult a hiányzó értékek kizárása után. A piacon töltött idő változót z-standardizáltam az együtthatók értelmezhetősége érdekében.

Két modellt hasonlítottam össze: hagyományos legkisebb négyzetek (OLS) regressziót és robusztus regressziót (MM-becslő), amely kevésbé érzékeny a kiugró értékekre.

Regressziós modellek létrehozása

# Regressziós modellek létrehozása - két modell: lm vs. lmrob (robusztus) 
mod_lm <- lm(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
  data = model_df)

mod_rob <- lmrob(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
  data = model_df)

# Regressziós modellek eredményei
summary(mod_lm)

## 
## Call:
## lm(formula = revision_number ~ atc_group + days_on_market_z + 
##     generic + orphan_medicine + accelerated_assessment, data = model_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.208  -4.405  -0.668   4.441  62.413 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 13.7840     0.7468  18.458  < 2e-16 ***
## atc_groupJ                   3.1040     0.9543   3.253  0.00118 ** 
## atc_groupL                   3.6263     0.8707   4.165 3.38e-05 ***
## atc_groupN                   2.2845     1.0407   2.195  0.02838 *  
## days_on_market_z             7.2030     0.3300  21.825  < 2e-16 ***
## genericTRUE                 -4.1782     0.8287  -5.042 5.46e-07 ***
## orphan_medicineTRUE         -1.0943     1.1333  -0.966  0.33449    
## accelerated_assessmentTRUE   3.6924     1.8587   1.986  0.04725 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.828 on 1018 degrees of freedom
## Multiple R-squared:  0.3808, Adjusted R-squared:  0.3765 
## F-statistic: 89.42 on 7 and 1018 DF,  p-value: < 2.2e-16

summary(mod_rob)

## 
## Call:
## lmrob(formula = revision_number ~ atc_group + days_on_market_z + generic + 
##     orphan_medicine + accelerated_assessment, data = model_df)
##  \--> method = "MM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.543  -3.935  -0.299   4.045  60.442 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 14.2598     0.4997  28.536  < 2e-16 ***
## atc_groupJ                   3.6720     0.8605   4.267 2.16e-05 ***
## atc_groupL                   2.9631     0.6137   4.828 1.59e-06 ***
## atc_groupN                   1.5387     0.6807   2.261  0.02399 *  
## days_on_market_z             8.8971     0.4217  21.099  < 2e-16 ***
## genericTRUE                 -3.5214     0.4465  -7.886 7.98e-15 ***
## orphan_medicineTRUE         -0.1621     0.5224  -0.310  0.75643    
## accelerated_assessmentTRUE   3.6881     1.1926   3.092  0.00204 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 6.065 
## Multiple R-squared:  0.6258, Adjusted R-squared:  0.6232 
## Convergence in 20 IRWLS iterations
## 
## Robustness weights: 
##  28 observations c(33,45,238,257,293,302,303,305,411,501,976,990,991,992,993,994,995,1014,1015,1018,1019,1020,1021,1022,1023,1024,1025,1026)
##   are outliers with |weight| = 0 ( < 9.7e-05); 
##  103 weights are ~= 1. The remaining 895 ones are summarized as
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001371 0.828100 0.953900 0.850500 0.986900 0.999000 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol          zero.tol 
##         1.000e-07         1.000e-10         1.000e-07         1.000e-10 
##       eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
##         9.747e-05         4.351e-12         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

Regressziós modellek összehasolítása

# Együtthatók összehasonlítása
coef_comparison <- bind_rows(
  tidy(mod_lm, conf.int = TRUE) %>% mutate(model = "OLS"),
  tidy(mod_rob, conf.int = TRUE) %>% mutate(model = "Robust"))

coef_comparison %>%
  select(model, term, estimate, std.error, conf.low, conf.high, p.value) %>%
  mutate(
    estimate = round(estimate, 3),
    std.error = round(std.error, 3),
    conf.low = round(conf.low, 3),
    conf.high = round(conf.high, 3),
    p.value = format_p(p.value))

## # A tibble: 16 × 7
##    model  term                     estimate std.error conf.low conf.high p.value
##    <chr>  <chr>                       <dbl>     <dbl>    <dbl>     <dbl> <chr>  
##  1 OLS    (Intercept)                13.8       0.747   12.3      15.2   < 0.001
##  2 OLS    atc_groupJ                  3.10      0.954    1.23      4.98  0.001  
##  3 OLS    atc_groupL                  3.63      0.871    1.92      5.34  < 0.001
##  4 OLS    atc_groupN                  2.28      1.04     0.242     4.33  0.028  
##  5 OLS    days_on_market_z            7.20      0.33     6.56      7.85  < 0.001
##  6 OLS    genericTRUE                -4.18      0.829   -5.80     -2.55  < 0.001
##  7 OLS    orphan_medicineTRUE        -1.09      1.13    -3.32      1.13  0.334  
##  8 OLS    accelerated_assessmentT…    3.69      1.86     0.045     7.34  0.047  
##  9 Robust (Intercept)                14.3       0.5     13.3      15.2   < 0.001
## 10 Robust atc_groupJ                  3.67      0.86     1.98      5.36  < 0.001
## 11 Robust atc_groupL                  2.96      0.614    1.76      4.17  < 0.001
## 12 Robust atc_groupN                  1.54      0.681    0.205     2.87  0.024  
## 13 Robust days_on_market_z            8.90      0.422    8.07      9.72  < 0.001
## 14 Robust genericTRUE                -3.52      0.447   -4.40     -2.65  < 0.001
## 15 Robust orphan_medicineTRUE        -0.162     0.522   -1.19      0.862 0.756  
## 16 Robust accelerated_assessmentT…    3.69      1.19     1.35      6.03  0.002

# Együtthatók vizualizációja
ggplot(coef_comparison %>% filter(term != "(Intercept)"), 
       aes(x = estimate, y = term, color = model)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), 
                 position = position_dodge(width = 0.5), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  labs(
    title = "Regressziós együtthatók összehasonlítása",
    subtitle = "OLS vs. Robusztus regresszió",
    x = "Becsült együttható (95% CI)",
    y = "Prediktor",
    color = "Modell"
  ) +
  theme_minimal()

Regressziós eredmények

OLS modell

Az OLS modell szignifikánsan illeszkedett az adatokra, F(7, 1018) = 89.42, p < .001, és a revíziószám varianciájának 38.1%-át magyarázta (R² = .381, R²adj = .377).

A piacon töltött idő volt a legerősebb prediktor: minden egységnyi növekedés a z-standardizált piacon töltött időben 7.20 egységgel növelte a revíziók számát, B = 7.20, 95% CI [6.56, 7.85], p < .001. Ez azt jelenti, hogy a régebben engedélyezett gyógyszerek több revízión estek át.

Az ATC csoport szignifikáns prediktornak bizonyult. A referenciakategóriához (A csoport) képest a J csoport (B = 3.10, p = .001), az L csoport (B = 3.63, p < .001) és az N csoport (B = 2.28, p = .028) egyaránt szignifikánsan magasabb revíziószámmal rendelkezett.

A generikus státusz szignifikáns negatív összefüggést mutatott a revíziók számával (B = -4.18, 95% CI [-5.80, -2.55], p < .001), ami arra utal, hogy a generikus gyógyszerek átlagosan 4.18-cal kevesebb revízión esnek át, mint az originális készítmények.

A gyorsított értékelés szignifikáns pozitív prediktora volt a revíziószámnak (B = 3.69, 95% CI [0.05, 7.34], p = .047). Az orphan gyógyszer státusz nem mutatott szignifikáns összefüggést (B = -1.09, p = .334).

Robusztus modell

A robusztus regresszió (MM-becslő) 28 megfigyelést azonosított kiugró értékként. A robusztus modell lényegesen magasabb magyarázóerővel rendelkezett (R² = .626, R²adj = .623), ami arra utal, hogy a kiugró értékek jelentősen torzították az OLS becsléseket.

A prediktorok mintázata hasonló maradt, de néhány fontos eltérés mutatkozott. A piacon töltött idő hatása erősebb volt a robusztus modellben (B = 8.90, 95% CI [8.07, 9.72], p < .001). Illetve a gyorsított értékelés hatása is robusztusabbá vált (B = 3.69, p = .002).

Modell összehasonlítás

A keresztvalidáció (80-20% train-test split) alapján a robusztus modell kissé jobb prediktív teljesítményt mutatott: RMSE = 11.5 vs. 11.7, MAE = 7.32 vs. 7.68, R²pred = .376 vs. .356.

A Cook-féle távolság alapján 51 megfigyelés (5.0%) minősült potenciális kiugró értéknek (Cook’s d > 4/n). A residuumok vizsgálata heteroszkedaszticitásra utalt, ami a robusztus módszerek alkalmazását indokolttá tette.

Regressziós modellek illeszkedése

# Modell illeszkedési mutatók
compare_performance(mod_lm, mod_rob, metrics = c("R2", "R2_adjusted", "RMSE", "AIC", "BIC"))

## # Comparison of Model Performance Indices
## 
## Name    | Model |    R2 |  RMSE |  AIC (weights) |  BIC (weights)
## -----------------------------------------------------------------
## mod_lm  |    lm | 0.381 | 9.789 | 7610.9 (>.999) | 7655.3 (>.999)
## mod_rob | lmrob | 0.626 | 9.956 |        (>.999) |        (>.999)

Regressziós modellek - Keresztvalidáció

# Regressziós modellek keresztvalidációja
set.seed(123)
idx <- sample(seq_len(nrow(model_df)), size = floor(0.8 * nrow(model_df)))
train <- model_df[idx, ]
test  <- model_df[-idx, ]

m1 <- lm(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment, 
  data = train)
m2 <- lmrob(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment, 
  data = train)

pred1 <- predict(m1, newdata = test)
pred2 <- predict(m2, newdata = test)

# Predikciós hibák
rmse <- function(y, yhat) sqrt(mean((y - yhat)^2))
mae  <- function(y, yhat) mean(abs(y - yhat))
r2_pred <- function(y, yhat) 1 - sum((y - yhat)^2) / sum((y - mean(y))^2)

prediction_comparison <- tibble(
  Modell = c("OLS", "Robust (lmrob)"),
  RMSE = round(c(rmse(test$revision_number, pred1), rmse(test$revision_number, pred2)), 3),
  MAE = round(c(mae(test$revision_number, pred1), mae(test$revision_number, pred2)), 3),
  R2_pred = round(c(r2_pred(test$revision_number, pred1), r2_pred(test$revision_number, pred2)), 3))
prediction_comparison

## # A tibble: 2 × 4
##   Modell          RMSE   MAE R2_pred
##   <chr>          <dbl> <dbl>   <dbl>
## 1 OLS             11.7  7.68   0.356
## 2 Robust (lmrob)  11.5  7.32   0.376

Regressziós modellek - Diagnosztika

# Regressziós modellek diagnosztikája (ábra)
check_model(mod_lm)

check_model(mod_rob)

Regressziós modellek - Kiugró értékek

# Cook's distance az OLS modellhez
model_df$cooks_d <- cooks.distance(mod_lm)
model_df$fitted <- fitted(mod_lm)
model_df$residuals <- residuals(mod_lm)

# Kiugró értékek azonosítása (Cook's d > 4/n)
n <- nrow(model_df)
outliers <- model_df %>%
  filter(cooks_d > 4/n)

cat("Kiugró értékek száma:", nrow(outliers), "\n")

## Kiugró értékek száma: 51

cat("Kiugró értékek aránya:", round(nrow(outliers)/n * 100, 1), "%\n")

## Kiugró értékek aránya: 5 %

# Vizualizáció
ggplot(model_df, aes(x = fitted, y = residuals)) +
  geom_point(aes(size = cooks_d, color = cooks_d > 4/n), alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("FALSE" = "gray50", "TRUE" = "red")) +
  labs(
    title = "Residuumok vs. Illesztett értékek",
    subtitle = "Piros pontok: kiugró értékek (Cook's d > 4/n)",
    x = "Illesztett értékek",
    y = "Residuumok",
    size = "Cook's d",
    color = "Kiugró"
  ) +
  theme_minimal()

Regressziós modellek - Összefoglaló táblázat

# Összefoglaló táblázat (modelsummary)

modelsummary(
  list("OLS" = mod_lm, "Robust" = mod_rob),
  estimate = "{estimate} ({std.error})",
  statistic = "p.value",
  stars = TRUE,
  title = "Regressziós modellek összehasonlítása",
  coef_rename = c(
    "atc_groupJ" = "ATC: J",
    "atc_groupL" = "ATC: L",
    "atc_groupN" = "ATC: N",
    "days_on_market_z" = "Piacon töltött idő (z-score)",
    "genericTRUE" = "Generikus",
    "orphan_medicineTRUE" = "Orphan gyógyszer",
    "accelerated_assessmentTRUE" = "Gyorsított eljárás"
  )
)

Regressziós modellek összehasonlítása
	OLS	Robust
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	13.784 (0.747)	14.260 (0.500)
	(<0.001)	(<0.001)
ATC: J	3.104 (0.954)	3.672 (0.860)
	(0.001)	(<0.001)
ATC: L	3.626 (0.871)	2.963 (0.614)
	(<0.001)	(<0.001)
ATC: N	2.285 (1.041)	1.539 (0.681)
	(0.028)	(0.024)
Piacon töltött idő (z-score)	7.203 (0.330)	8.897 (0.422)
	(<0.001)	(<0.001)
Generikus	-4.178 (0.829)	-3.521 (0.447)
	(<0.001)	(<0.001)
Orphan gyógyszer	-1.094 (1.133)	-0.162 (0.522)
	(0.334)	(0.756)
Gyorsított eljárás	3.692 (1.859)	3.688 (1.193)
	(0.047)	(0.002)
Num.Obs.	1026	1026
R2	0.381	0.626
R2 Adj.	0.377	0.623
AIC	7610.9
BIC	7655.3
Log.Lik.	-3796.426
F	89.422
RMSE	9.79	9.96

Összegzés és konklúzió

Jelen elemzés az Európai Gyógyszerügynökség gyógyszer-engedélyezési adatbázisát vizsgálta négy fő kutatási kérdés mentén.

Engedélyezési trendek. A gyógyszer-engedélyezések száma jelentős fluktuációt mutatott 1995 és 2023 között, kiugró növekedéssel a COVID-19 pandémia időszakában. Az L csoport (daganatellenes szerek) dominanciája különösen az elmúlt évtizedben vált hangsúlyossá.

ATC csoportok közötti különbségek. Szignifikáns különbségek mutatkoztak a revíziók számában (H(3) = 9.39, p = .025) és a piacon töltött időben (H(3) = 48.2, p < .001) a négy leggyakoribb ATC csoport között. Az L csoport gyógyszerei szignifikánsan rövidebb ideje vannak a piacon, ami feltételezhetően az onkológiai területen tapasztalható gyorsabb innovációval magyarázható.

Revíziószámot befolyásoló tényezők. A regressziós elemzés alapján a piacon töltött idő, az ATC csoport, a generikus státusz és a gyorsított értékelés szignifikáns prediktorai a revíziók számának. A robusztus regresszió magasabb magyarázóerőt mutatott (R² = .63 vs. .38), ami a kiugró értékek jelentős torzító hatására utal.

Hálózati kapcsolatok. A központi idegrendszeri gyógyszerek hálózatelemzése feltárta a terápiás területek és hatóanyagok közötti komplex kapcsolatrendszert. A Parkinson-kór, a Szkizofrénia és az Epilepszia központi pozíciót foglal el a hálózatban, ami ezeknek a betegségterületeknek a terápiás diverzitását tükrözi.

Limitációk

Az elemzés korlátai közé tartozik, hogy az adatbázis csak az EMA által centralizált eljárásban engedélyezett gyógyszereket tartalmazza, így az eredmények nem általánosíthatók a nemzeti szinten engedélyezett készítményekre. A revíziók száma mint kimeneti változó korlátozottan tükrözi a gyógyszerek tényleges biztonságossági és hatékonysági profilját. A hálózatelemzés keresztmetszeti jellegű, és nem tárja fel az időbeli változásokat a kapcsolatrendszerben.

Jövőbeli irányok

További kutatások vizsgálhatnák a revíziók tartalmát és típusát, a forgalomból történő kivonások prediktorait, valamint az engedélyezési folyamat időtartamának meghatározó tényezőit.

Statisztikai programozás - ELTE PhD - Final project

Girasek Hunor (RW8V2Q)

2025-12-15

Bevezetés

Adatok előkészítése

Adatok feldolgozása

Minta

Gyógyszer-engedélyezési trendek

ATC-rendszer

Az ATC-rendszer szerinti megoszlás

Statisztikai elemzések - TOP4 ATC

Leíró statisztika

A normalitás vizsgálata

Csoportok összehasonlítása

Revíziók száma

Piacon töltött idő

Hálózatelemzés

Adatok előkészítése

Szűrés és hálózat építése

Hálózat vizualizáció

Hálózatemezés mutatók

Centralitás elemzés

Regressziós modellek

Regressziós modellek létrehozása

Regressziós modellek összehasolítása

Regressziós eredmények

OLS modell

Robusztus modell

Modell összehasonlítás

Regressziós modellek illeszkedése

Regressziós modellek - Keresztvalidáció

Regressziós modellek - Diagnosztika

Regressziós modellek - Kiugró értékek

Regressziós modellek - Összefoglaló táblázat

Összegzés és konklúzió

Limitációk

Jövőbeli irányok